%% input signal
global M;
global N; %bits for quantifying the coefficients
N = 10; % coefficients quantization
M = 16; % input signal quantization
Fs = 16000;
Fin = 500;                   
T = 1/Fs;                     % Sample time
NSTART = 1000;                % Length of signal
NFFT = 8192;
L = NSTART + NFFT;
Ncycles = max(primes(NFFT*Fin/Fs));
Fin = Ncycles * Fs/ NFFT;
t = (0:L-1)*T;   
x =  1*sin(2*pi*Fin*t); 
XFFT = fft(x(NSTART+1:L),NFFT);

xq = (round((2^(M-1)-1)*x)) / 2^(M-1);
xqw = xq(NSTART+1:L) .* blackmanharris(NFFT)';
% xqw = xq;
XFFTQ = fft(xqw,NFFT);

%% filter information
Wp = 3400/8000; 
Ws = 4600/8000;
Rp = 0.3;
Rs = 60;

%% ideal
[n,Wp] = ellipord(Wp,Ws,Rp,Rs);
[num,den] = ellip(n,Rp,Rs,Wp);
y = filter(num,den,xq);
y = y(NSTART+1:L) .* blackmanharris(NFFT)';
YFFT = fft(y, NFFT);

%% transfer function form
[h,w] = freqz(num,den,NFFT,Fs);
[z,p,k] = tf2zp(num,den);

zq = round((2^(N-1)-1)*z) / 2^(N-1);
pq = round((2^(N-1)-1)*p) / 2^(N-1);
kq = round((2^(N-1)-1)*k) / 2^(N-1);

[numz,denz] = zp2tf(zq,pq,kq);
numq = round((2^(N-1)-1)*numz) / 2^(N-1);
denq = round((2^(N-1)-1)*denz) / 2^(N-1);
[hq,wq] = freqz(numq,denq,NFFT,Fs);
ytf = filter(numq,denq,xq);
ytf = ytf(NSTART+1:L) .* blackmanharris(NFFT)';
YFFTtf = fft(ytf,NFFT);

%% sos without scaling 
[sos,g] = zp2sos(zq,pq,kq);
gq = round((2^(N-1)-1)*g) / 2^(N-1);
ys = gq*sosfilter(sos(3,1:3), sos(3,4:6),sosfilter(sos(2,1:3),sos(2,4:6),sosfilter(sos(1,1:3),sos(1,4:6),xq)));
ys = ys(NSTART+1:L) .* blackmanharris(NFFT)';
YFFTs = fft(ys, NFFT);

%% infinity-norm scaling with up-ordering
[sosinf,ginf] = zp2sos(zq,pq,kq,'up','inf'); 
gqinf = round((2^(N-1)-1)*ginf) / 2^(N-1);
yinf = gqinf*sosfilter(sosinf(3,1:3), sosinf(3,4:6),sosfilter(sosinf(2,1:3),sosinf(2,4:6),sosfilter(sosinf(1,1:3),sosinf(1,4:6),xq)));
yinf = yinf(NSTART+1:L) .* blackmanharris(NFFT)';
YFFTinf = fft(yinf, NFFT);

%% 2-norm scaling with up-ordering
[sos2,g2] = zp2sos(zq,pq,kq,'up','two');  
gq2 = round((2^(N-1)-1)*g2) / 2^(N-1);
y2 = gq2*sosfilter(sos2(3,1:3), sos2(3,4:6),sosfilter(sos2(2,1:3),sos2(2,4:6),sosfilter(sos2(1,1:3),sos2(1,4:6),xq)));
y2 = y2(NSTART+1:L) .* blackmanharris(NFFT)';
YFFT2 = fft(y2, NFFT);

%% SNR
widthl = 5;
widthr = 5;
fmin = 0;
fmax = Fs/2;
snrx    = snr(XFFT,   NFFT,Fs,widthl,widthr,fmin,fmax)
snrxq   = snr(XFFTQ,  NFFT,Fs,widthl,widthr,fmin,fmax)
snry    = snr(YFFT,   NFFT,Fs,widthl,widthr,fmin,fmax)
snrytf  = snr(YFFTtf, NFFT,Fs,widthl,widthr,fmin,fmax)
snrsos  = snr(YFFTs,  NFFT,Fs,widthl,widthr,fmin,fmax)
snrinf  = snr(YFFTinf,NFFT,Fs,widthl,widthr,fmin,fmax)
snr2    = snr(YFFT2,  NFFT,Fs,widthl,widthr,fmin,fmax)

%% plot
f = Fs*linspace(0,1,NFFT);

subplot(4,2,1);
semilogx(f,20*log10(abs(XFFT)));
grid on;

subplot(4,2,2);
semilogx(f,20*log10(abs(XFFTQ)));
grid on;

subplot(4,2,3);
semilogx(f,20*log10(abs(YFFT)));
grid on;

subplot(4,2,4);
semilogx(f,20*log10(abs(YFFTtf)));
grid on;

subplot(4,2,5);
semilogx(f,20*log10(abs(YFFTs)));
grid on;

subplot(4,2,6);
semilogx(f,20*log10(abs(YFFTinf)));
grid on;

subplot(4,2,7);
semilogx(f,20*log10(abs(YFFT2)));
grid on;



